/*
 #
 #  Files       : Vecteur.h
 #                ( C++ header file )
 #
 #  Description : The SmallMatrix Library
 #                ( http://code.google.com/p/smallmatrix )
 #
 #  Copyright   : Olivier Juan
 #                ( http://www.mas.ecp.fr/vision/Personnel/juan/ )
 #
 #  License     : CeCILL-C
 #                ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #  This software is governed by the CeCILL-C license under French law and
 #  abiding by the rules of distribution of free software.  You can  use,
 #  modify and or redistribute the software under the terms of the CeCILL-C
 #  license as circulated by CEA, CNRS and INRIA at the following URL
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL-C license and that you accept its terms.
 #
*/

#ifndef _VECTEUR_H
#define _VECTEUR_H
#include <SmallMatrix/SmallMatrix.h>

NAMESPACE_BEGIN(SmallMatrix)

template <int N> class Vecteur : public Matrix_Base<Vecteur<N>,N,N,1>
{
public:
	Vecteur(): Matrix_Base<Vecteur<N>,N,N,1>(){};
	Vecteur(const reel tab[N]):Matrix_Base<Vecteur<N>,N,N,1>(tab){}
	template <typename TT,int SSIZE> Vecteur(const Matrix_Base<TT,SSIZE,N,1>& copy):Matrix_Base<Vecteur<N>,N,N,1>(copy){}
	Vecteur(const reel &v):Matrix_Base<Vecteur<N>,N,N,1>(v){} // Pour tout N
	Vecteur(const reel &x,const reel &y);					  // N=2 seulement
	Vecteur(const reel &x,const reel &y,const reel &z);		  // N=3
	inline reel& operator()(const int& idx,const int& null=0);
	inline const reel& operator()(const int& idx,const int& null=0)const;
	template <typename T1,int SIZE1,int MN,typename T2,int SIZE2> inline Vecteur<N>& prod(const Matrix_Base<T1,SIZE1,N,MN>& M1,const Matrix_Base<T2,SIZE2,MN,1>& M2);
	template <typename T1,int SIZE1,typename T2,int SIZE2> inline Vecteur<N>& add(const Matrix_Base<T1,SIZE1,N,1>& M1,const Matrix_Base<T2,SIZE2,N,1>& M2);
	template <typename T1,int SIZE1,typename T2,int SIZE2> inline Vecteur<N>& subst(const Matrix_Base<T1,SIZE1,N,1>& M1,const Matrix_Base<T2,SIZE2,N,1>& M2);
	template <typename T1,int SIZE1> inline Vecteur<N>& add(const Matrix_Base<T1,SIZE1,N,1>& M1);
	template <typename T1,int SIZE1> inline Vecteur<N>& subst(const Matrix_Base<T1,SIZE1,N,1>& M1);
	template <typename T1,int SIZE1> inline Vecteur<N>& mult(const Matrix_Base<T1,SIZE1,N,1>& M1);
	template <typename T1,int SIZE1> inline Vecteur<N>& div(const Matrix_Base<T1,SIZE1,N,1>& M1);
	inline reel dot(const Vecteur<N>& Two)const;
	template <int P> inline Matrix<N,P> outer_product(const Vecteur<P>& Two)const;
	inline SquareMatrix<N> outer_product(const Vecteur<N>& Two)const;
	inline Vecteur<3> operator^(const Vecteur<3>&v) const;
	inline bool operator<(const Vecteur<N>& test)const;
	inline Vecteur<N>& normalize();
};

template <> 
inline Vecteur<2>::Vecteur(const reel &x,const reel &y):Matrix_Base<Vecteur<2>,2,2,1>(){
	(*this)(0)=x;(*this)(1)=y;
}

template <> 
inline Vecteur<3>::Vecteur(const reel &x,const reel &y,const reel &z):Matrix_Base<Vecteur<3>,3,3,1>(){
	(*this)(0)=x;(*this)(1)=y;(*this)(2)=z;
}

#pragma warning( disable : 4100)
//! Read/Write Accessor
template <int N> inline reel& Vecteur<N>::operator ()(const int& idx,const int& null)
{
	assert((idx>=0)&&(idx<N)&&(null==0));
	Vecteur<N>& ref=*this;
	return ref[idx];
}
//! Read Only Accessor
template <int N> inline const reel& Vecteur<N>::operator ()(const int& idx,const int& null)const
{
	assert((idx>=0)&&(idx<N)&&(null==0));
	const Vecteur<N>& ref=*this;
	return ref[idx];
}
#pragma warning( default : 4100)

//! Matrix-Matrix Product in which the result is a Vecteur
template <int N> template <typename T1,int SIZE1,int MN,typename T2,int SIZE2> inline Vecteur<N>& Vecteur<N>::prod(const Matrix_Base<T1,SIZE1,N,MN>& M1,const Matrix_Base<T2,SIZE2,MN,1>& M2) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		for (int k=0;k<MN;k++)
			ref(i,0)+=M1(i,k)*M2(k,0);
	return ref;
}
//! Matrix-Matrix Sum in which the result is a Vecteur
template <int N> template <typename T1,int SIZE1,typename T2,int SIZE2> inline Vecteur<N>& Vecteur<N>::add(const Matrix_Base<T1,SIZE1,N,1>& M1,const Matrix_Base<T2,SIZE2,N,1>& M2) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		ref(i,0)=M1(i,0)+M2(i,0);
	return ref;
}
//! Matrix-Matrix Substract in which the result is a Vecteur
template <int N> template <typename T1,int SIZE1,typename T2,int SIZE2> inline Vecteur<N>& Vecteur<N>::subst(const Matrix_Base<T1,SIZE1,N,1>& M1,const Matrix_Base<T2,SIZE2,N,1>& M2) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		ref(i,0)=M1(i,0)-M2(i,0);
	return ref;
}
//! Matrix-Matrix Sum in which the result is a Vecteur
template <int N> template <typename T1,int SIZE1> inline Vecteur<N>& Vecteur<N>::add(const Matrix_Base<T1,SIZE1,N,1>& M1) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		ref(i,0)+=M1(i,0);
	return ref;
}
//! Matrix-Matrix Substract in which the result is a Vecteur
template <int N> template <typename T1,int SIZE1> inline Vecteur<N>& Vecteur<N>::subst(const Matrix_Base<T1,SIZE1,N,1>& M1) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		ref(i,0)-=M1(i,0);
	return ref;
}
template <int N> template <typename T1,int SIZE1> inline Vecteur<N>& Vecteur<N>::mult(const Matrix_Base<T1,SIZE1,N,1>& M1) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++)
		ref(i,0)*=M1(i,0);
	return ref;
}
template <int N> template <typename T1,int SIZE1> inline Vecteur<N>& Vecteur<N>::div(const Matrix_Base<T1,SIZE1,N,1>& M1) {
	Vecteur<N>& ref=*this;
	for (int i=0;i<N;i++) {
		assert(M1(i,0)!=0);
		ref(i,0)/=M1(i,0);
	}
	return ref;
}
//! Inner Product with the Vecteur \a Two
template <int N> inline reel Vecteur<N>::dot(const Vecteur<N> &Two)const
{
	const Vecteur<N>& ref=*this;
	return ref[N-1]*Two[N-1]+(*((const Vecteur<N-1>*)(this))).dot(*((const Vecteur<N-1>*)(&Two)));
}
template <> inline reel Vecteur<1>::dot(const Vecteur<1> &Two)const
{
	const Vecteur<1>& ref=*this;
	return ref[0]*Two[0];
}
//! Outer Product (*this) * Two'
template <int N> template <int P> inline Matrix<N,P> Vecteur<N>::outer_product(const Vecteur<P> &Two)const
{
	const Vecteur<N>& ref=*this;
	Matrix<N,P> A;
	for (int i=0; i<N; i++) for (int j=0; j<P; j++) A(i,j)=ref(i)*Two(j);
	return A;
}
template <int N> inline SquareMatrix<N> Vecteur<N>::outer_product(const Vecteur<N> &Two)const
{
	const Vecteur<N>& ref=*this;
	SquareMatrix<N> A;
	for (int i=0; i<N; i++) for (int j=0; j<N; j++) A(i,j)=ref(i)*Two(j);
	return A;
}
template <int N> inline bool Vecteur<N>::operator<(const Vecteur<N>& test)const
{
	const Vecteur<N>& ref=*this;
	bool res=true;
	int i=0;
	while((res)&&(i<N))
	{
		res=(ref[i]==test[i]);
		i++;
	}
	i--;
	return (res?false:(ref[i]<test[i]));
}
template <> 
inline Vecteur<3> Vecteur<3>::operator^(const Vecteur<3>&v) const {
		Vecteur<3> res;
		res(0)=(*this)(1)*v(2)-(*this)(2)*v(1);
		res(1)=-(*this)(0)*v(2)+(*this)(2)*v(0);
		res(2)=(*this)(0)*v(1)-(*this)(1)*v(0);
		return res;
}
template <int N> inline Vecteur<N>& Vecteur<N>::normalize() 
{
	return ((*this) /= Matrix_Base<Vecteur<N>,N,N,1>::norm());
}

NAMESPACE_END

#endif
